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Abstract 

Cosmic rays impacting on the atmosphere cause particle-showers. 
Several descriptions exist for the evolution of the shower size along 
the atmospheric depth. The well known functions for shower profiles, 
Greisen, Gaisser-Hillas and 'Gaussian in Age', are intimately connected 
in that they all are approximate solutions of versions of the Rossi 
and Greisen diffusion equations. The mathematical connection will 
be demonstrated by means of two simple models for the longitudinal 
electromagnetic shower profile. Both models can be regarded either 
as a generalization of the Heitler model or as a simplification of the 
diffusion model of Rossi and Greisen. These models are far closer to 
reality than the Heitler model, while they are not as close to reality as 
the model of Rossi and Greisen. Therefore, they will be referred to as 
intermediate models. For each intermediate model the evolution of the 
shower is governed by either a single differential equation or a single 
integro-differential equation. The approximate solution of the differ- 
ential equation is a Gaisser-Hillas function and can be adjusted such 
that it almost matches the Greisen profile. The approximate solution 
of the integro-differential equation is a 'Gaussian in Age' function. The 
corresponding profile is, after suitable adjustment, in excellent agree- 
ment with the Greisen profile. The analysis also leads to an alternative 
functional form for the age parameter. 
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1 Introduction 



The longitudinal development of electromagnetic showers can be described 
by a system of diffusion equations. They can be solved by means of functional 
transforms and a second order saddlepoint approximation [TJ [2l [3]. The 
solution of Rossi and Greisen can be elaborated to what is known as the 
Greisen function [3]. For this some further approximations had to be made. 
Because of the inaccuracies involved in these approximations one may ask if 
a satisfying trial function can also be obtained from a less accurate approach 
[5j. This is valid as long as the inaccuracy in the model does not disturb too 
much the essential shape of the profile. Fortunately, as we will see, it does 
not. Deviations in height and width can be easily adjusted for. In this way 
we obtain a simple route to the construction of trial functions for shower 
profiles. 

In this paper we will restrict ourselves to the longitudinal development of 
the electromagnetic cascade governed by the three elementary processes of 
pair production, Bremsstrahlung and ionization losses. The corresponding 
shower consists of three particles: electrons, positrons and photons. A simple 
model for the longitudinal development of the electromagnetic cascade is the 
Heitler model [6j. More recently, the Heitler model has been applied for 
the development of the hadronic portion in the initial stages of extensive 
air showers [7]. According to the Heitler model each particle will split after- 
travelling the same distance into two particles of half the energy of the parent 
particle. This distance is the splitting length d = A r ln2, where A r is the 
radiation length in the atmosphere: A r ~ 37.1 g/cm 2 and d ~ 25.7 g/cm 2 . 
We start at atmospheric depth X = 0. After the first collision there will 
be 2 particles at atmosferic depth X = d. The subsequent collisions lead 
to 4 particles at X = 2d, 8 particles at X = 3d and so on. At atmospheric 
depth nd there will be 2 n particles in the shower. At this depth the energy of 
each particle is Eq/2 71 , where E$ is the energy of the primary particle. This 
cascade continues until the energy of the particles falls below the critical 
value E c = 84 MeV, the energy at which the ionization loss is equal to 
the collisional energy loss. The shower then stops, according to the Heitler 
model, when n > n c , where n c = In (Eq/E c ) / In 2. 

2 An intermediate shower model 

In reality particles do not travel equal distances before they split. To model 
this we cut the atmosphere into slices AX of equal atmospheric thickness. 
Travelling through a slice, a step from now on, each particle with a certain 
energy has a chance p to split into 2 particles of half that energy, and a 
chance q = 1 — p to continue as a single particle with the original energy. 
This is a generalization of the Heitler model. After n steps the energies E(k) 



2 



the particles can posses are 



E(k) = ±E , (1) 



where only the levels with < k < n can be occupied. After n steps the 
expected number of particles in level k is N (k, n) . The conservation of energy 
then requires 

£iiV(M) = l- (2) 

k=0 

Since a fraction q remains in the energy level and a fraction 2pN enter from 
the higher energy level the expectation values N(k, n + 1) are related to the 
N(k,n) as follows 

N(k,n + 1) = qN(k,n) + 2pN(k- l,n). (3) 

This relation obeys energy conservation: 

n -, n+1 -, 

E^«) = 1 ^2 ¥ N(k,n + l) = l. (4) 

k=0 k=0 

As can be verified by induction, the solution of eqs. (|3]) and |4|) is 



n 



This is the same as having 2 k particles of energy E$/2 k with binomial prob- 
ability 



n 



P(k,n)=^ k jp*q n -«. (6) 

That is, the splittings are binomially distributed, while k splittings lead to 
2 k particles. The total expected number of particles after n steps is 

N ^ n ) = E (f\ (2p) k <l n ~ k = (2p + q) n = (1 + P) n - (7) 



k=0 



Thus N(n) initially grows exponentially. The expected fraction of particles 
at energy level k after n steps then is 

N(k,n) _ fn\ (2p) k (l- P r- k 

This is actually a binomial distribution with redefined probability: 

P^J^ => f(k,n)=( n \p) k (l-pr- k . (9) 
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As in each splitting the energy is conserved, E n = Eq, the average energy 
per particle after n steps is 

. , E(n) E Q , . 

en = = T, — (10) 

v ' N{n) (l+p) n v ; 

For the average number of steps it takes before a particle splits we consider 
the probability wi for a particle to survive I steps, but not I + 1 steps: 

wi = q l p. (11) 

These probabilities are also properly normalized: 

oo oo 

E™i=pEv 1 = t^ = 1 - (12) 

1=0 1=0 q 

For the average number of steps a particle survives without splitting we then 
obtain 

OO l OO i -i 

//^ Y^ ; d V- i 1 PI 1 1 i / r Q^ 

<l>=2^lwi= pq— 2^q = PQ^-z. = 71 ^2 = - = - - 1. 13 

T^O dq i=1) dql-q {l-q) 2 p p 

Since (< I > +1) times AX is the splitting length d and since n times AX is 
the actual atmospheric depth X, we obtain for small values of p the following 
relation: 

AY" Y 

d=(<l>+l)AX = = — . (14) 

p np 

For p = 1 the original Heitler model is retained: X = nd. Our aim is to 
consider the situation in the limit p — )■ 0. 

3 Continuum limit 

In the continuum limit, p — > and n — > oo satisfying (|14p . the difference eq. 
Q turns into a differential equation which can be written either as 

9N &n ^ = qN<Kkl ^ + 2pN ( k ~ !' n ) ~ N ( k > n ) = 2 P N ( k - 1j «) - pN(k, n) 

(15) 

Since the binomial distribution in the limit p — > equals the Poisson distri- 
bution, one can expect solutions of the form 2 k times a Poisson distribution. 
Indeed the equations (|15p and (fl~6j) are solved by 

N(k,n) = ^-(2np) k e~ np (17) 
fc! 
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and 

1 / ? X \ k x 

^• x) = l(-) ^ < 18 > 

respectively. These solutions satisfy the boundary condition N(0, 0) = 1. 
Notice that the number of particles in energy level E(k) = E^j2 k reaches 
its maximum at depth X max = kd. Since d = A r In 2 this is X max = 
\ r ln(Eo/E). The expression for N(k,n) or N(k,X) is the same as hav- 
ing 2 k particles of energy Eo/2 k with Poisson distributed probability: 

P(k,n) = y(np) k e- n P, (19) 

or 

p <^>4(f) Vf - < 2o > 

In other words, the splittings are Poisson distributed with average length d, 
while k splittings lead to 2 k particles. We will refer to (fT8|) as the Poisson- 
related distribution. 

Without absorption the total expected number of particles at atmospheric 
depth X, thus after n steps, is 



^) = El^)Vf ,21, 



which in practice becomes 



(¥)*«■*=«*■ (22) 

That is, without absorption the total number of particles grows as e x l d . 



4 Absorption 

The shower model presented above is applicable to the early stages of the 
shower, when all particles have relatively large energy. At later stages a 
particle will be absorbed or scattered out of the shower when its energy 
drops below the critical value E c ~ 84 MeV. Then the energy of the shower 
is no longer conserved. The net growth of the number of particles will slow 
down. After reaching a maximum the number of particles will decrease and 
finally the shower will fade out as far as the particles have not reached the 
surface of the earth yet. Since the energy of a particle after k splittings is 
Eq ■ 2~ k , the particle is taken out of the shower when k > [*n c ], where the 
critical parameter n c is given by 

n c =^m. (23) 
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In the discrete model we take the ceiling m = \n c ] as the stopping value 
since a particle with energy Eq = 2 nc ■ E c may split m times before the 
energy is below E c . As a consequence the average total number of particles 
develops in time as 

n 

N(n) = Y^0{m- k)2 k P(k,n), (24) 

A:=0 

where 6(x) is the Heaviside step function 

*(*) = { n v X -n ( 25 ) 
v ' 1 0, if x < 0. y J 

The previous result for N(n) holds as long as n < m: 

N(n< m] X) = ±±(^p) k ''e- x ' d . (26) 
However, for n > m the summation is limited: 

N(n>m-,X) = pl(^ k e-^ (27) 

Since n = goes to infinity in the limit p — > 0, the latter equation can in 
practise be used at all depths: 

W = tl( 2 4) k ^ X/d - (28) 



k =o kl V d J 



It can also be written as 



A m 

N m (X) = - ■ e x ' d ■ f(X; k + 1, d/2) (29) 

k=0 

or as 

_ x/d T(m + l,2X/d) 

Nm[x) ~ 6 r(m + i) ' (30) 

where / is the Gamma distribution 

X k-i . -x/e 

m^)= et . m ■ pi) 

where T(m+ 1, 2X/d) is the incomplete Gamma function and where T is the 
Gamma function: T(m + 1) = m! if m is an integer. By taking the derivative 
of expression (|28p with respect to X we obtain 
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as well as 



<iX <i dm\\d 

From eqs. (|16|) and (|32|) we see that the total number of particles satisfies 
the same differential equation as the number of particles in each energy level. 
This suggests to take a single term of the Poisson series as an approximation 
for the total number of particles. Numerical evaluation of the expression (|28j) 
demonstrates that the total number of particles reaches a maximum A max 
at atmospheric depth 

TTL^ d 

X max « -r « (m - l)d. (34) 

m + 1 

Therefore the forelast term of the series (1281) seems to be the best choice: 



N m {JC]nJ^(^>T\- x l* i (35) 



r(m) \ d 

where A is a normalization constant. This approximation actually is a 
Gamma distribution. It is known that a truncated Poisson series can be 
accurately approximated by a Gamma distribution. Another motivation for 
the Gamma distribution lies in the fact that it satisfies the same differen- 
tial equation as the truncated Poisson series. We will make further remark 
about our motivation for the Gamma distribution at the end of section 8. 
As desired (|35|) reaches its maximum at A max = (m — l)d. From eqn. (|33|) 
we find: 

iV max - mT{m) y d J e , (3oj 

while for the approximation (|35|) this is 



A /2X 



m— 1 

AW « =^r ( ) e-*— / d . (37) 

r(m) V a 



Comparison of the latter two eqs. gives 



A«4(l-i + ofi 5 ]J. CIS) 



m V ^ 2 



By means of this expression for ^4 and by means of the Stirling approxima- 

/ 1 \ m— 1 

tion, r(m) ~ ylis^m — 1 ( ^^j 1- J , the function (|35p takes the form 



y™^l \( m - l )dJ 



where 

1 



4:r 



^)«!-t: + o( j). (40) 
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Notice that g is close to unity. As an example g ~ 0.97 for m = 25. As known 
the agreement between a Gamma distribution and a truncated Poisson series 
can be further improved by means of a parameter Xq as follows: 

vr / m _ l \(m - l)o - x y 

Accurate fits are obtained with Xq ~ —d/2. The larger the negative values 
for Xq, the larger the width of the profile, the fuuhm for instance. Writing 
= as iV max , (m — l)d as X maa; and d as A, the expression (|4ip reads 



m- 2 



AT m (X) ~ ^max T ^ e s . (42) 

We clearly recognize it as the Gaisser-Hillas function [8]. For the moment 
we will restrict ourselves to the case Xq = — ^d. We will adjust the width of 
the profile afterwards. The approximation of a truncated Poisson series by 
a Gamma distribution is pure mathematics. The parameter Xq therefore is 
just a mathematical parameter. This (and the fact that its value is prefer- 
ably negative) supports the opinion that it should not be given a physical 
interpretation as being the point of first interaction [9) 110). 

The function m = [~n c ] approximately goes as m ~ n c + \. The substi- 
tution of m ~ n c + i in the expression ([41 I) leads to 

N m « 5 (n c ) AZl (l+^\ n \n^,2-X /d , (43) 



To visualize its accuracy both the numerical summation (|28p and the function 
(|43p are plotted in figure 1. The expression (|43|) leads to the same profile as 
the following exression 

jv m (x) « 4=4= (4i ) e " c ~ x/d > ( 44 ) 



7r ,/n7 Vn,,d 



except that it is translated over a distance d/2. In the next section we will 
make a comparison with the Greisen function. Since we are interested in 
the main characteristics of the shape of the profile and not in the differ- 
ence caused by a small translation, we will conveniently use expression 
hereafter. 



5 The ratio of the numbers of particles 

Usually one works with an expression for the number of charged particles in 
the shower instead of the total number of particles. To this end we will dis- 
tinguish the electrons and positrons from the photons. An electron/positron 
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Vertical is the number of particles and horizontal is the atmospheric depth. 



can split into an electron/positron and a photon, while a photon splits into 
an electron and a positron. It doesn't change the present model, while we 
can now look at the following system of difference equations: 

N e ±(n + 1) = N e ±(n) + 2piV 7 (n), 
iV 7 (ra + 1) = P N e ±{n) + (1 - p)N 7 (n). (45) 

Here N e ± is the number of electrons and positrons and iV 7 is the number of 
photons. For the ratio p = N e ±/N^ we find 

p(n + l)iV 7 (n + 1) = p(n)N,y(n) + 2pN y (n), 

N 7 (n + 1) =pp(n)JV 7 (n) + (l-p)JV 7 (n). (46) 

The elimination of iV 7 leads to 

p(n) + 2p 

P(« + I)=i — 7 r- (47) 

1 — p + pp(n) 

The ratio /) asymptotically approaches a limit value R, which satisfies R 2 — 
R — 2 = 0. The latter equation is solved by the stable stationary point 
R = 2. From iV e± = 2iV 7 we obtain iV 7 = ±iV and iV e± = |iV. To obtain 
the number of charged particles we therefore have to multiply the expression 
(|44p by |. So, under the assumption that the shower develops with the ratio 
in its equilibrium according to the intermediate model the number of charged 
particles is given by 

4 2 nc (eX\ nc x 



A re±m = 9 („ c )__^_ 1 (48) 
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6 Expressions for the shower profile 



The shower profile according to the expression (j48J) turns out to be twice 
too high, twice too narrow compared to, for instance, the Greisen profile. A 
deviation from the Greisen profile could be expected since in reality there 
will also be other than fifty-fifty splittings with different probabilities for 
different ratios. Moreover, for Bremsstrahlung these probabilities differ from 
the ones for pair creation. If this is taken into account, one arrives at the 
system of diffusion equations of Rossi and Greisen jT] . We will return to it in 
the next section. Here we will modify the eq. (|48p such that the height and 
width of the corresponding profile is in agreement with the Greisen profile. 
We halve the height by dividing (|48p by 2. We double the width of the profile 
by multiplying the powers in (|48p by | In 2. The result is 

2 2 nc feX\l ncla2 



N e± (X)=g{n c )- T —- — -e~i3-. (49) 

\n c d; 



This still is a Gaiser-Hillas type of function as can be verified by substituting 

3d 
2 In 2 

also be written as 



X m ax = n cd, A = an< ^ -^0 = ^ ex P ress i° n (02]). The eqn. (|49l) can 



N e± (t) = g(n c )-^^- ■ (^-) .e-f, (50) 

where t = j- is the atmospheric depth in units of radiation length. Its 
maximum value 

2 2™ c 

-ZV e ±,max = g(n c )—=—= (51) 
3V 71 " y/n c 



occurs at depth t max = n c ln2. 

The well-known Greisen approximation formula reads [2]: 

N e± (t) = ^.e^-I^), (52) 
Dc 



where s = f+ 3 ^ is the age-parameter and y c = In (Eq/ E c ) = n c ln2. Com- 
pletely in terms of t and n c this is 

0.31 (\ 2n c ln2\§* , . . 

Its maximum value 

31 2™ c 2™ c 
iV e±imax = -= • — = w 0.37 • -— (54) 
V In 2 V^c V^c 

also occurs at depth i max = n c \a.2. For a visual comparison both the present 
Gaisser-Hillas profile and the Greisen profile are plotted in the figure 2. It 
shows that both profiles practically coincide. From the comparison of (|5ip 
with (|54p we conclude that the value 0.31 in the Greisen function approxi- 
mately equals the semi-theoretical value 2 ^? . 
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Figure 2: Longitudinal shower profiles according to the present Gaisser- 
Hillas function (solid) and the Greisen function (dashed) for a 2 • 10 15 eV 
shower. Vertical is the number of electrons and positrons and horizontal the 
atmospheric depth (in units of radiation length). 



7 The connection with Rossi and Greisen 

In this section we will show the connection between the preceeding model and 
the cosmic-ray theory of Rossi and Greisen PQ. For our purposes it suffices to 
consider the situation under what is known as the 'approximation A'. Under 
this approximation the complete screening cross sections of radiation and 
pair creation processes are used. Other processes like Compton scattering 
are neglected and ionization loss is solely used as a stopping criterion [TJ [3] . 
The diffusion equations for the differential distributions n e ± and n 7 then 
read 

f°° 1 f E' E\ r E 1 fW\ 

j E n e± {E',t)-<p Q [— —)dE'-J Q ne± (E,t)-<p [—)dW, (55) 

-i;^ E ^ (D « - f «^<>> (I) «* 

(56) 

In these equations (fo {i^j "jjdVK i s the differential probability per radiation 
length for an electron (or a positron) with energy between E and E + dE 
to split of a photon with energy between W and W + dW. Similarly, 
V>o (w 7 ) * s ^ ne differential probability per radiation length for a photon 

with energy between W and W + dW to produce a pair with the electron 
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energy between E and E + dE. They are given by [T] 



V-oO-t) = t^o(l/«) = \u 2 - \u + 1 + 26(n 2 - u). (57) 

The last term on the rhs is of minor importance since b is relatively small, 
b ~ 0.013. After a suitable change of variables the equations <j55[) and (|56p 
take the form 



<9t Jo \u ) J U J(j \ 1 — D ' y 1 

-n e± (E,t) / <^o(v)df, (58) 
Jo 

an 7( w, t) _ /-i ng± ^ ^) du _ n ^ t) £ Mu)du m 



dt 

By means of separation of the variables energy and depth one obtains the 'el- 
ementary' solutions: n e ±(E,t) = aE~^ s+l ^>e xt and n 7 (W, t) = bW~^ s+1 ^e xt . 
Upon substitution a quadratic equation for A is obtained. One solution, 
A2(s), corresponds to a quick adaptation to the equilibrium ratio between 
n e ± and n 7 . The other solution, Ai(s), describes the main development of 
the shower, see section 27 of reference [T|. The connection of the intermedi- 
ate model with the model of Rossi and Greisen follows when delta-functions 
are substituted for the probabilities. The eqns (|55p and (|56p then take the 
form 



dn e ±(E,t) 
dt 



2 / n 7 (W, t) a5(W - 2E)dW + 
Je 



roo rE 

/ n e ± (E', t) a5{E' - 2E)dE' - n e ±(E, t) / a5{2W - E)dW, (60) 

JE JO 



dn^(W,t) 



roo rW 

/ n e ±(E,t)a5(E - 2W)dE -nJW,t) / ad(2E - W)dE, 
Jw Jo 

(61) 

where a = j— g since the differential probabilities in the diffusion equations 
are per radiation length A r , while the splitting probability is a delta- function 



dt 



per splitting distance d = A r ln2. The system (|60p and ()6ip reduces to 
dn e ±(E,t) 



Ot 



2an 1 (2E, t) + an e ±{2E, t) - an e ±(E, t) (62) 



dn 7 (W,t) 

— - ( ( / I , ± \ _ M . I I — II I I - 



an e± (2W,t)-anJW,t). (63) 



In the hypothetical one-particle model there also is no difference between 
photons, electrons or positrons. Given the 1 : 2 ratio for the number of 
photons and electrons/positrons in our model, we take n e ±(E,t) = ^n(E,t) 
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and n.y(W,t) = ^n(E,t). We can also set W equal to E. Then the system 
of equations further reduce to 

2dn(E,t) 2 . , 2 . , 2 . 

3 m = ^an(2E, t) + -an(2E, t) - -an(E, t) (64) 

ldn(E,t) 2 . . 1 , 

3 dt = ^n(2E, t) - -an(E, t). (65) 

Obviously, these equations are identical. As it should, the one-particle model 
is governed by a single differential equation: 

Q = 2an(2E, t) - an(E, t). (66) 

Substituting a = we S e ^ 

dn(E,t) 2 . . 1 . , 

_L^i = _„( 2£ , ( )-_„( E , t ), (67) 

or 

^^l = 2 -n(2E,X)- 1 -n(E,X). (68) 
From = E • 2~ fc it follows that E(k - 1) = 2E(fe). Hence, 

*» I| ^,W t -l),A-)-l, Wt ),I). (69) 
This is the same differential equation as in (|16p . 



8 Another intermediate model 

In the previous model we considered splittings with the energy equally di- 
vided between the decay particles. We obtain a little more accurate model 
by allowing other than fifty-fifty splittings, although all with equal probabil- 
ity. That is, we take the differential probabilities for the different splittings 
equal to a constant. It turns out that the maximum shower size occurs at 
the desired depth if the constant is taken equal to 2. So, for the following 
analysis we will restrict ourselves to the case ipo(u) = y>o{v) = 2. Then the 
system of equations (|58p and (|59|) take the form 

dn e ±(E,t) f 1 (E \ 2 , f 1 ( E \ 2 , 
-2 I n J (—,t)-du+l n e ±[- -,*)- -dv 



dt Jo \u J u Jo V 1 — v J 1 — v 

-n e ±(E,t) f 2dv, (70) 
J o 



dn y (W,t) 
dt 



r 1 (W \ 2 r 1 

' n e ±\—,tJ-dv-n~,(W,t)J 2du. 



(71) 
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After a suitable change of variables this is: 
dn e ±(E,t) 



dt 

dn 7 (W,t 



/■l f E \ 2 f E \ 2 

2/ n 7 (— ,i ) -du+ / n e ± I — ,t ) ~du-2n e ±(E,t), (72) 
Jo ' \u J u Jo \u J u 

f 1 (W \ 2 

/ n e ± — ,t ) -du-2nJW,t). (73) 
Jo \ u J u 



Also in this model we do not distinguish between photons and charged par- 



tides. By setting n e ±(E,t) = 2n^(W, t) = %n(E,t) and W = E, we obtain 



2dn(E,t) 4 f 1 (E \ 1 , 4 Z" 1 \ 1 4 .„ , 

' n — ,t -dn + -/ n — ,t ) - - -n(E,t) (74) 



(75) 



3 9f 3 7o \u J u 3 7o \u J u 

1 dn(E,t) 



A f 1 f E \ I 2 .„ , 
3 Jo \u J u 3 



Obviously, these equations are identical. Also this one-particle model is 
governed by a single diffusion equation: 



dn(E,t) 
dt 



4 j 1 n(—,t] -du-2n(E,t). (76) 
Jo \u J u 



This equation allows for elementary solutions in which the variables E and 
t are separated. To be specific, solutions of the type 

n(E,t) = A.(^-Y +1 -e xt , (77) 



. E . 

with s —1, do satisfy the differential equation (|76[) if 

A( S ) = 7^-2. (78) 

Notice that if s does depend on t, we would have the additional require- 
ment y + X't = 0, where y = 1ii(Eq/E) and where the prime stands for the 
derivation with respect to s. This requirement is precisely the saddle point 
condition (|87p as we will see soon. Now we will construct the solution in 
an analogous manner as in the paper of Rossi and Greisen. To this end we 
consider the Mellin integral 

M n (s,t)= E s n(E,t)dE (79) 
Jo 

and its inverse transformation 

1 rc+ioo 

n{E,t) = — E- s ^M n (s,t)ds. (80) 

Z7TI J c—ioo 



14 



Multiplying both sides of eqn. (|76p by E s and integrating with respect to 
energy from to oo, we obtain: 



oo 



at Jo Jo \u J u 

P oo 

2 / E s n(E,t)dE. (81) 
J o 



Since n(E/u,t) = u^ s+1 ^ n(E , t) , the latter is reduced to 



-/ E s n(E, t)dE = — — / (S) s n(E,t)d(E)- 
at Jo s + 1 Jo 

/•oo 

2 / E s n(E,t)dE. (82) 
Jo 



Hence, 



^M„(s,t) = A(s)M„(s,t), (83) 
with A(s) as given by (|78p . The solution of this equation is 

M n (s,t) = M n (s,0)-e x ^ t , (84) 

where 

roo roc 

M n {s,0)= E s n(E,0)dE = E S 5(E — E )dE = Eq, (85) 
Jo Jo 

Next we apply the inverse Mellin transformation: 

rc+ioo 



1 rc+too 

n(E,t) = — E~ s ^M n (s,t)ds 

Z7VI Jc—ioo 
1 rc-\-ioo 

= — ; / E- s - l E s e x{s)t ds 
QiTri J c— zoo 

1 1 [t (Eq\ x(s)t 



2m' Eq Jc-ioo V E 
1 1 rc+ioo 

~ E~o~2~rd 



ds 

e y.( S+ l)+X( S )t d ^ (gg) 



where y = \u{Eq/E) = A; In 2 is the lethargy. The real constant c must be 
in the strip of analyticity, which is the positive halfplane. We already saw 
in section 3 that the distribution in energy level E reaches a maximum at 
tmax = V- At this depth the exponent y-(s+l)+X(s)t is equal to t max -(s+l) + 
X(s)t max . For s along the real axis it has a minimum at s given by A'(s) = — 1, 
That is, for s = 1 and thus A = 0. Since an analytic function satisfies the 
Cauchy-Riemann equations and thus the Laplace equation, the exponential 
term should have a maximum at the point s, along directions perpendicular 
to the real axis. Although one usually does not know in advance the relation 
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between y and t max , the foregoing makes it clear that we can require the 
integrand to have a saddle point at the point s defined by 

y + \'(s)t = 0, (87) 

where the prime stands for differentiation with respect to s. A second order 
Taylor series of the exponent of the integrand around the point s then yields 

y ■ (s + 1) + X(s)t ^ys + y + \{s)t + -X"(s)t(s - s) 2 . (88) 
Taking the integration path through the saddlepoint, we obtain 

Eq 2iri Js-ioo 
With the change of variables, s = s + ix, this is 

n {E,t) = ±-^- e y°+y+mt f°° e -^"(^ 2 dx. (90) 

Eq 2-7T 7_oo 

Evaluating the Gaussian integral, we obtain 

1 e ys+y+Ks)t 

Next we require the solutions to reach a maximum. The function e ys+y+> *^ t 
reaches its maximum at depth t max given by 

[y + X'(s)t max ] W + A(s) = 0. (92) 

Because of the saddlepoint relation (|87p the latter implies A(s) = 0. From 

(|78p it is inferred that the maximum occurs at s = 1. Since A'(l) = — 1 it 

follows from (|87j) that t max = y or X max = \ r \n(Eo/E) as desired. From 
(|57|) and d7HD it follows that 

\"(s)t = y^ (93) 

\{s)t = - 2t (94) 

and 

S - = 2^-l. (95) 

Hence, 

n(E,t) = — e - 7 ==. (96) 

Eq V 27r • V 
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For the integral distribution we can also make use of the Mellin transform: 

poo 

M N (s-l,t)= E S ~ 1 N(W > E,t)dE «-> 
Jo 

1 rc+ioo 

N(W > E,t) = — / E~ s M N (s - 1, t)ds. (97) 

-^7T2 Jc— zoo 



As can be verified by means of partial integration and the property ^- f(x)dx 

— f(y), there holds the following relation between the Mellin transforms of 
the integral and differential distribution: 

M N (s-l,t) = -M n (s,t), (98) 
s 

where N = N(W > E,t) = JJ° n(E ,t)dE' . From this relation we obtain 

1 rc+ioo 1 

N{W > E,t) = — / -E- s M n (s,t)ds. (99) 

Z7TI J c—ioo S 

Substitution of the expression (fM|) for M n leads to 

— [ C+i0 °- (^-Ye^ds = — [ C+i °° e -^s+ys + X(s)t dS} (10Q) 
2-7TZ Jc-ioo s \ E J 2-iri J c -ioo 

where again 

A( s ) = -L--2. (101) 

From here we can proceed in a similar manner as for the differential distri- 
bution. The integrand has a saddle point at the point s defined by 

^- + y + X'(s)t = 0. (102) 

s 

A second order Taylor series of the exponent of the integrand around the 
point s then yields 

- Ins + y • s + X(s)t « - Ins + y ■ s + X(s)t + -r(s,i)(s - s) 2 , (103) 
where 

r(l,t) = (^ + A"(l)t). (104) 

Taking the integration path through the saddlepoint, we obtain to second 
order 

N(W >E,t) = — : / ' e -^s+ys+mt+r(s,t)±(s-sf ds _ ^ 1Q5 ) 

27T2 J s _2oo 

With the change of variables, s = s + ix, this is 

N(W > E,t) = -L e -^+ys+X(s)t / e - 5 r(S,t)* dx _ ( 1Q6 ) 
27T ./-co 
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Evaluating the integral, we obtain 

p — In s+y-s+\(s)t 

N(W>E,t) = . (107) 

V27rr(s,t) 

Next we require the solutions to reach a maximum. The function e~ lns +f 
reaches its maximum at depth t max given by 



-- + y + X'(s)t r 
s 



ds\ 

tit) 



+ A(s) = 0. (108) 



Because of the saddlepoint relation (I102|) the latter implies A(s) = 0. From 
(jlOip it is inferred that the maximum occurs at s = 1. It follows from (jlOip . 
(fT02|) and (fTUD that _ 

r(s,t)^y^, (109) 
A(s)t « - 2t (110) 



and 



Hence, 



^> E ,*) = ^f.-^. (112) 

First we will consider the situation without the factor \J . At the end of 

this section we will show that the influence of this factor can be neglected. 

For E = E c and thus y = y c = n c ln2 we obtain for the age parameter, 
leaving the bar, 

8 = 2* — I 1 (113) 

V n c \u2 v ' 

and for the shower size 



g 4Vin c ln2-2i 



N(t) = N(W > E c , t) = = . (114) 

2 • 2 nc V 2ir ■ n c In 2 



For the charged part this is 



e 4Vin c ln2-2i 



Ne±(t) = ——=_=. (115) 

3 • 2 nc ^2ir • n c m2 



The latter can also be written as 

0.16 • 2 r 



N e ±(t) = — e -^mctt^ncinij i / n6 ) 
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By means of the age parameter (|113p it can also be written as 

Ne±(t) = 2^_J!! e -^ c i n 2( S -i )2 (117) 



n c 

We clearly recognize it as a Gaussian in Age. Notice that for the present age 
parameter s = — 1 if t = 0. Although the corresponding profile does have the 
right shape with a maximum at the right position, t max = ra c ln2, it differs 
in height and width from the Greisen profile. To obtain the same height we 
simply replace 0.16 by 0.37. To obtain the same width we replace the factor 
^ in the exponent by ^. Then the eqs. (|116p and (QT7J read 

Ne± (t) = Q - 37 " 2 "V f(*-2vCT2+n e ln2) (118) 



and 



Ne± ( t ) = 0-37-2" e e _i ncln2( ,_ 1)2 ( | lf)) 



respectively. From equation (|5ip it is inferred that the numerical value 0.37 
in the expression for the shower size is almost equal to the semi-theoretical 
value: -^m- So, we can also write the Gaussian in Age profile as 

N e± (t) = J—^-e-^ 2 ^' 2 . (120) 
3V 71 " V^c 



The latter Gaussian in Age profile has standard deviation a = , = , by 

y ZiTlc 111 ^ 

means of which the Gaussian in Age profile can also be written as 

2 2 nc 1 1\2 
N e ± (t) = . = • e-i(— ) . 121 

In figure 3 both the present Gaussian in Age profile and the Greisen 
profile are plotted. We see the profiles match nicely. For practical purposes 
the Gaussian in Age profile can be generalized to a three parameter trial 
function 

N e± (t) = N max ■ e --(^-v^) . (122) 

The parameter w determines the width of the profile; its value will be close 
to |. As for the Gaiser-Hillas function, it can be generalized further to a four 
parameter function by means of the shift t — y t — £q and t max — ^ t max — ^o- 



Finally we will consider the situation where the factor J y | is not neglected. 
Then the expression (|118p should be modified accordingly. That is, (I118p 
should be muliplied by a factor f nc j n2 j 

Ne±{t) = ^L£l f^V\ e -U^Vi^ + n c ln2) (123) 
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2.0 xlO 6 




-> t 



Figure 3: Longitudinal shower profiles according to the present 'Gaussian 
in Age' function (solid) and the Greisen function (dashed) for a 2 • 10 15 eV 
shower. Vertical is the number of electrons and positrons and horizontal the 
atmospheric depth (in units of radiation length). 



2.0 xlO 6 



1.5 xlO 6 - 

+i 

Z 1.0 xlO 6 - 

r ; 

500000 - 

- 




10 20 30 40 

-> t 

Figure 4: Longitudinal shower profiles according to the expression (|123|) 
(solid) and the Greisen function (dashed) for a 2 • 10 15 eV shower. Vertical is 
the number of electrons and positrons and horizontal the atmospheric depth 
(in units of radiation length). 
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In figure 4 both the profile according to <j!23|) and the Greisen profile are 
plotted. The influence of the factor is small and practically equal to a shift 
over a distance \d. As mentioned before, a small translation is not of interest 
for the present comparison of profiles and will therefore be ignored. 

Although redundant a similar analysis as in this section can in principle 
be applied to the eq. (|69[) as well. Such an analysis leads exactly to the 
solution (|18|) for the number of particles in the discrete energy levels and 
approximately to the Gamma distribution (|35j) . In fact this was our main 
motivation to take the Gamma distribution for the total number of particles 
in section 4. 



9 Summary and conclusions 

To day both the Gaisser-Hillas function and the Gaussian in Age function 
are used as trial functions for the reconstruction of longitdinal shower profiles 
[T2] . There are even trial functions composed of two halves of Gaussian 
in Age functions [13J. All the these function usually contain parameters 
which are not independent of each other [14]. Recently new less dependent 
parameters for the Gaisser-Hillas and the Gaussian in Age functions were 
constructed |14| I15j. By introducing the parameter [i = y c = n c ln2, the 
expressions for the shower profile the espression for the shower profile can be 
written by means of a single parameter. The Greisen function, the Gaisser- 
Hillas function and the Gaussian in Age function then respectively read 



0.31 (\ 2fi 



2' 



0.31 ft\h 



N e± (t) = ^-[-)' -eM*, (125) 
and 

AT e± (t) = 2^1 • e - W§ v^-f*. (126) 

In this minimal form all three profiles practically coincide and reach a max- 
imum • at depth t max = [i. If desired one can replace the numerical 

constant 0.31 by its semi-theoretical analogon 2v J^ . It is readily admit- 
ted that we crudely neglected small translations and ignored the differences 
between individual shower profiles in order to obtain the single parameter ex- 
pressions. Of course, to take account for small translations and other details 
of observed or simulated individual shower profiles additional parameters are 
unavoidable. 

The fact that the three profile functions practically coincide gives rise to 
the idea that there must be a mathematical connection or common origin. 
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We found the connection by solving the Rossi and Greisen equations for 
simplified cross sections. The similarity of the three profiles also leads to 
the conclusion that the shape of the shower profile is rather independent of 
the type of cross section. Instead we conclude that the characteristic shape 
of the shower profile is governed by the statistics of the splittings. Only 
splittings with substantial different probability distributions, such as in the 
Heitler model, will lead to a substantial different profile. The functional form 
of the cross sections mainly influence the height and the width of the profile. 

Because of its mathematical convenience the Gaiser-Hillas function is 
sometimes favoured over the Greisen function |16| . The Gaussian in Age 
function is mathematically convenient as well. The Greisen function is less 
convenient. The points of inflection, for instance, are for the Gaiser-Hillas 

function the roots fi ± of the quadratic equation t 2 — 2 fit + /j, 2 — |/x = 0. 
For the Gaussian in Age function the points of inflection are two of the roots 
of the cubic equation t 3 — 2/jt 2 + (fi 2 — — ^/j = 0. The degree of these 
polynomials expresses the hierarchy in the complexity of the corresponding 
expressions for the shower profiles. For the Greisen function the points of 
inflection can not be derived analytically. 

In trial functions for longitudinal shower profiles the Greisen age param- 
eter 

a = - (127) 

t 4- 2/ 

is commonly used. However, one should be careful with regarding it as a 
universal age parameter since the functional form of the age parameter is 
model dependent [17J. With the use of the parameter fj, the Gaussian in Age 
function reads 

N e± (t) = °-^-^.e-^-V 2 , (128) 
while the present analysis suggests 

- 1 (129) 

as the natural age parameter. It seems worthwhile to investigate the alter- 
native age parameter on its practical use. This subject is currently under 
research. 




22 



Acknowledgements 

I am grateful to the reviewers for their valuable suggestions and improve- 
ments. I wish to thank Prof. J.W. van Holten and Prof. B. van Eijk for their 
encouragement and support. I also wish to thank Nikhef for its hospitality. 
The work is supported by a grant from FOM (Foundation for Fundamental 
Research on Matter). 

References 



[I] B. Rossi, K. Greisen, Rev. Mod. Phys. 13, 240 (1941). 

[2] J.F. Carlson, J.R. Oppenheimer, Phys. Rev. 51, 220 (1937). 

[3] J. Nishimura, Handbuch der Physik, XLVI/2, 1 (1967). 

[4] K. Greisen, Prog. Cosmic Ray Phys. 3, 1 (1956). 

[5] R.W. Schiel, J.P Ralston, Phys. Rev. D 75, 016005 (2007). 

[6] W. Heitler. The Quantum theory of Radiation, Oxford Univ. 

Press (1954). 

[7] J. Matthews, Astropart. Phys. 22, 387 (2005). 

[8] T. Gaisser, A.M. Hillas, Proc. 15th ICRC, Plovdiv, Bulgaria 

8, 353 (1977). 

[9] C. Song, Astropart. Phys. 22, 151 (2004). 

[10] P. Sommers, Compt. Rendus Phys. 5, 463 (2004). 

[II] M. Unger et al., Nucl. Instr. Meth. Phys Res. A588, 433 
(2008). 

[12] T. Abu-Zayyad et al, Astropart. Phys. 16 1 (2001). 

[13] M. Giller et al., J. Phys. G: Nucl. Part. Phys. 31, 947 (2005). 

[14] J.A.J. Matthews et al, J. Phys. G: Nucl. Part. Phys. 37, 

025202 (2010). 

[15] S. Andringa et al, Astropart. Phys. 34, 360 (2011). 

[16] J. Linsley, 19th ICRC, NASA Conf. Publ. 7, 167 (1985). 

[17] P. Lipari, Nucl. Phys. B (Proc. Suppl.) 196, 309 (2009). 



23 



